A numerical study of two-photon ionization of helium using the Pyprop framework 
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Few-photon induced breakup of helium is studied using a newly developed ab initio numerical framework 
for solving the six-dimensional time-dependent Schrodinger equation. We present details of the method and 
calculate (generalized) cross sections for the process of two-photon nonsequential (direct) double ionization at 
photon energies ranging from 39.4 to 54.4 eV, a process that has been very much debated in recent years and is 
not yet fully understood. In particular, we have studied the convergence property of the total cross section in the 
vicinity of the upper threshold (~54.4 eV), versus the pulse duration of the applied laser field. We find that the 
cross section exhibits an increasing trend near the threshold, as has also been observed by others, and show that 
this rise cannot solely be attributed to an unintended inclusion of the sequential two-photon double ionization 
process, caused by the bandwidth of the applied field. 

PACS numbers: 32.80.Rm, 32.80.Fb, 42.50.Hz 



I. INTRODUCTION 

The study of how light interacts with matter has occupied 
physicists for a long time. Of particular fundamental inter- 
est is single and multiple ionization of atoms and molecules 
by photon impact, with subsequent ejection of one or a multi- 
ple number of electrons. In this respect, single-photon multi- 
ple ionization is special, as exchange of energy between the 
involved electrons is a prerequisite for the process to take 
place. The investigation of such correlated dynamical pro- 
cesses poses many unique challenges to experiment and the- 
ory. A prime example of this is the one-photon double ion- 
ization of helium, which has been the subject of intense study 
since the pioneering work of Byron and Joachain As a 
matter of fact, it was only recently that a complete agreement 
between theoretical calculations and accurate measurements 
with synchrotron radiation was established, for the value of 
the (generalized) cross section for the direct (nonsequential) 
double ionization process 10-111] • 

The problem of two-photon double ionization of two- 
electron atomic systems presents additional difficulties. First, 
the separation between sequential and nonsequential double 
ionization often becomes a nontrivial problem [60, in situa- 
tions where both processes are energetically accessible |01- 
Here, 'sequential' ionization usually refers to the process 
where the electrons are emitted one after the other by sub- 
sequent absorption of a photon each, and where the second 
electron has time to relax to a stationary ionic state before it is 
emitted. Thus, energy exchange between the two electrons is 
not strictly required. In contrast, 'nonsequential' double ion- 
ization depends upon exchange of energy between the outgo- 
ing electrons, and as such it represents a clear departure from 
an independent-particle picture. As mentioned above, after 
more than 40 years of investigation, the case of double elec- 
tron ionization by single photon absorption is now considered 
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to be a well understood process J^, but the related problem 
of two-photon direct (nonsequential) double ionization in he- 
lium, at photon energies ranging from 39.4 to 54.4 cV, is still 
being debated ||5L I9l422ll . In particular, the separation of two- 
photon single ionization, where the remaining electron is left 
in an (excited) ionic state, from two-photon double ionization, 
has turned out to be a subtle theoretical problem, as the role of 
electron correlations in the final states is not yet fully under- 
stood lUt]. Moreover, moving beyond the single-photon ion- 
ization regime is extremely challenging from the experimen- 
tal point of view, due to the weakness of the signals. Although 
the total cross section for the two-photon nonsequential dou- 
ble ionization of helium has been the subject of experiments, 
employing state of the art high-order harmonic ll23i l24ll and 
free-electron laser (FEL) radiation ifzsll . the experimental re- 
sults remain inconclusive ll26ll . 

In this paper, we revisit the problem of two-photon nonse- 
quential double ionization in helium. We present details on a 
recently developed B-spIine based numerical framework for 
solving the two-electron time-dependent Schrodinger equa- 
tion. The method is built on the more general Pyprop frame- 
work IztIi and was recently used to study the role of electron 
correlations in the stabilization of helium in intense extreme- 
ultraviolet (XUV) laser fields [28]. The two-electron mod- 
ule we have implemented is designed to utilize massively 
parallel supercomputers to perform accurate large-scale time- 
dependent simulations, and we here use it to make a contribu- 
tion to the ongoing discussion on two-photon double ioniza- 
tion. 

For the total cross section we obtain values that are in close 
agreement with some recently reported results ifTsi l2Tll . In 
contrast, our results are in clear disagreement with results 
from calculations where correlation effects in the final scat- 
tering states have been included to some extent, using differ- 
ent approximative methods. In our approach, which is similar 
to the method adapted by Feist et al. fistj, and others, this 
problem is circumvented by letting the ionized wave packet 
propagate for some time after the end of the laser pulse, so 
that most of the wave packet eventually reaches the asymp- 
totic (Coulomb) region ll29ll . The wave function is then ana- 
lyzed by means of projections onto a set of uncorrelated eigen- 
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states, i.e.. Slater determinants constructed from one-electron 
orbitals. 

With our present method we are able to consider relatively 
long pulses (up to about 10 fs in total pulse duration). This 
puts us in a position to study the convergence property of 
the total cross section in the vicinity of the upper threshold 
(~ 54.4 cV), as a function of pulse length. It has been re- 
ported that the cross section exhibits an apparent sharp rise 
near the threshold (H [M [H 111 S ■ This rise stands some- 
what in contrast to the results obtained in lowest (nonvanish- 
ing) order perturbation theory (LOPT) |6]. Our calculations 
with longer pulses indicate that the cross section indeed ap- 
pears to grow as it approaches the threshold (up to 51.7 eV), 
but we cannot rule out the possibility that the cross section 
reaches a maximum at some point in the immediate neigh- 
borhood of (or on) the threshold. Based on our results we 
are tempted to conclude that the increase of the cross section 
around 54.4 eV is not solely due to an unintended inclusion 
of the sequential process laH^. If this is correct, the very 
interesting question remains: what is the underlying physical 
process causing the unexpected behavior? 

The rest of this paper is organized into two main sections. 
The first outlines some theoretical background, and then goes 
on to discuss the various aspects of the numerical method, 
ending with a discussion of convergence issues. In the final 
section, we present our calculations of the total cross section, 
and pay particular attention to its behavior near the sequential 
threshold. 

Atomic units (Ti, nig and e replaced by 1) are used in the 
following exposition, except where otherwise noted. 



II. THEORY AND NUMERICAL APPROACH 



B-splines Bill have long been popular in atomic and 
molecular physics computations, due to their ability to accu- 
rately represent atomic eigenstates (see ll32ll and references 
therein). For time-dependent calculations, integration of the 
Schrodinger equation directly in the B-spline representation 
is very efficient for one-electron systems, due to the sparse 
and structured matrices that arise f3y\. For two-electron sys- 
tems, however, the matrix structure becomes more compli- 
cated, and the matrix sizes are also much larger Therefore, 
in time-dependent approaches for two-electron systems using 
B-spline discretization, one-electron orbitals or atomic eigen- 
states have been used to construct a matrix representation of 
the field interactions. Both these approached are useful and 
accurate, but do not easily scale to very large basis sizes, be- 
cause the basis functions are global, resulting in dense matri- 
ces. Operations involving such matrices are difficult to paral- 
lelize efficiently, which eventually becomes necessary as the 
basis size is increased. 

In this section, we present some details of a recently devel- 
oped B-spline based numerical approach to the two-electron 
time-dependent Schrodinger equation, where the time integra- 
tion is performed directly in the B-spline basis. A Python/C-H- 
implementation of the method have been created, which uses 
the Pyprop framework ll27ll . 



A. Theoretical background 

We consider the two-electron time-dependent Schrodinger 
equation (TDSE), 
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Employing the semi-classical approximation, the light-atom 
interaction Hamiltonian can be cast into the form. 
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where Hfi represents the interaction with the external field. 
The laser-atom interaction is modeled in the dipole approx- 
imation using the velocity gauge formulation, which, when 
linearly polarized along the z axis, can be written as 



H 



(3) 



Here is the vector potential defining the external field. 
The corresponding electric field is given by Ez = —dAz/dt. 
For the temporal form of the laser interaction, a sine-squared 
carrier-envelope was chosen, i.e.. 



Az{t) = Aosin-^ { — ] cos(wt), 
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(4) 



where Aq — Eq /uj, Eq being the electric field amplitude, u 
is the (central) laser frequency, and T defines the (total) pulse 
duration. We have also considered pulses with a Gaussian 
envelope. 



Az{t) = Aoexp 
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where tq is the full-width at half-maximum (FWHM) pulse 
duration, and T — 2to defines the (chosen) total pulse dura- 
tion. 



B. Discretization 

Turning the continuous TDSE into a set of coupled ordinary 
differential equations is achieved by representing the wave 
function, in spherical coordinates, as a product of radial B- 
splines and Coupled Spherical Harmonics, 
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(6) 

Here, k = {L, M, li, I2} is a combined index for the angular 
indices, and the Coupled Spherical Harmonics are given in 
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terms of Spherical Harmonics as ll34l] 

- Y.^hhmM - m\LM)Yl^{n,)Y,l'--\n^). (7) 

m 

For the special case of z-polarization, the problem reduces 
effectively to five dimensions, as M is conserved during the 
laser-atom interaction. Since we are studying ionization from 
the ground state {M = manifold), only 7\f = is in- 
cluded in the calculations. The B-spline basis functions de- 
pend upon several parameters, and determining the optimal 
values of these are not trivial (see Bachau et al. ll32il for a 
discussion). Throughout this paper, we have used order 5 B- 
splines, and an exponential distribution of breakpoints near 
the origin, with linear spacing further out, providing accu- 
rate representation of both bound and continuum states. Zero 
boundary conditions are imposed by removing the first and 
last B-spline for each radial direction. 

Except for the electron-electron interaction, all the terms 
in Eq. (O are one-electron operators, and are straightforward 
to discretize. When calculating matrix elements of these, the 
radial integrals are performed with Gauss-Legendre quadra- 
ture |j32[], while the angular integrals are handled analytically. 
The electron-electron interaction is first expanded in a trun- 
cated multipole series, 

i:7^-E i: ^^^^V.n{^^)Yu.[n,i (8) 

1=0 m=-l > 

and each of the terms are handled in a similar manner to the 
one-electron operators. Finally, taking into account the non- 
orthogonality of the B-spline basis, j drBi{r)Bj{r) = Sij, 
the TDSE can be written in matrix form 

js|^c(t)=H(t)c(i). (9) 
ot 

The total overlap matrix S is a Kronecker product of one- 
electron overlap matrices for each angular momentum com- 
ponent, 

S = Ifc ® Si (g> Sa- (10) 

In contrast to the one-electron case, the Hamiltonian matrix H 
does not have a simple banded structure, although it is quite 
sparse. An overall banded structure remains, but is now inter- 
laced with bands of zeros. 

At present, we make use of multiprocessor systems by dis- 
tributing the wavefunction across several processor in the an- 
gular rank, such that all the time-independent radial matrix 
elements are processor-local. This of course restricts the num- 
ber of processors we may use to the total number of angular 
momentum elements. We are currently implementing the op- 
tion of distributing one or both radial ranks, which would al- 
low us to use more processors, and thus an even larger radial 
basis. 



C. Time integration 

It is common for numerical integrations schemes to be 
based on the first-order exponential approximation to the 
propagator, 

T{U,tf) = exp[-iS-^nAt] +0{At'^), (11) 

which is quite accurate for reasonably small time steps At = 
tf — ti. The matrix exponential may be calculated efficiently 
by the popular Arnoldi or Lanczos methods Ha a. 
However, instabilities prevent us from using this approach in 
the present case; instead, we use the implicit Cayley-Hamilton 
form of the propagation operator, 

(^S + c(i + At) =(^S- c{t). (12) 

The above linear system of equations is typically too large 
to be solved by direct methods, but very sparse. It is there- 
fore solved using the iterative Generalized Minimum Residual 
method (GMRES) SHI. Similarly to the Arnoldi method, 
GMRES uses a Krylov subspace, constructed from succes- 
sive multiplications of (S + ^H) on c(i) in each time step. 
A least-square problem is solved in the subspace spanned by 
these vectors, and a solution of the equation with a mini- 
mum residual is obtained. With this method, the error in the 
computed solutions (residual) is controlled by the size of the 
Krylov subspace, which can be increased automatically, thus 
always ensuring a high precision solution. 

As is typical for discretized partial differential equations, 
and in particular those obtained with a B-spline basis, the 
system is quite stiff, and a good preconditioner is essential 
for the convergence of GMRES. Let A = (S + ^H) and 
b = (S — ^H)c. The preconditioner M is then constructed 
to make the linear system 

M-iAc = M-ib (13) 

easier to solve than the original system. This can be achieved 
by letting M be an approximation of A. The preconditioner 
used in this paper is a block type preconditioner, where each 
block consists of the complete radial Hamiltonian for a given 
coupled spherical harmonic (k is the angular index), 

M(i^j^k),{i' ,j' ,k')^(s{ij,k),(i' ,j' ,k) + —H(,^j^k),(i',j',k)jSk,k'- 

(14) 

This block diagonal matrix is distributed across processors in 
such a manner that the elements in the wave function corre- 
sponding to one block are all local to one processor When 
solving linear systems involving M, each block can be solved 
separately and thus there is no communication between dif- 
ferent processors. Furthermore, as the exact solution of M 
is not required, we employ the incomplete LU (ILU) factor- 
ization 1 38], M w LU, as provided by the IFPACK toolkit 
available in the Trilinos |l39] library. A similar preconditioner, 
using SuperLU 14^ to solve M exactly, was also tested, but 
found to be less efficient for the given system. 
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D. Extracting physical information 

Although excitation of the neutral atom is negligible for the 
intensities and frequencies we will consider in this paper, cal- 
culation of a subset of the eigenstates of the helium atom is 
nevertheless useful in many cases. The implicitly restarted 
Arnoldi method (IRAM) |41] may be used for this purpose. 
For reasons similar to those prompting the use of a precondi- 
tioner above, IRAM will converge slowly for interior eigen- 
values. However, shifted inverse iterations can be used to ac- 
celerate the convergence for eigenvalues near a given shift a. 



(H-aS)'^c, 



1 



£■„ - cr 



c„S. 



(15) 



IRAM requires the multiplication of the operator matrix on a 
vector in order to operate. For inverse iterations this corre- 
sponds to solving the linear system in Eq. (fTSl l, for the matrix 
B = (H — crS). Note the similarity between A and B, the dif- 
ference being only the scalars in front of H and S. The linear 
solver used for the propagation can therefore also be used to 
find eigenvalues with IRAM. Alternatively, a preconditioned 
Davidson method can be used. We found that when our basis 
grows sufficiently large, the Block-Davidson approach |42ll . 
as implemented in the Trilinos package Anasazi, performed 
favorably compared to the shift-invert Arnoldi method. In 
most cases, either of these methods can be used to rapidly 
obtain eigenpairs in the vicinity of any given shift value. 

In order to extract double ionization probabilities from the 
wavefunction, we must project onto some set of states which 
span this space. The exact double continuum states are hard to 
find, as they require solving a scattering problem for the full 
two-particle system. An approximation using non-correlated 
states, obtained by solving a set of one-electron radial eigen- 
value problems, is used instead. 



19^ 
29^ 
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The double ionization continuum is represented by a product 
of He+ (Z — 2) states, which, when expanded in the B-spline 
basis, we will denote b„j n^.fc. These states are not orthogo- 
nal to the bound states of the atomic system, and consequently 
the projection of the final wave function on the atomic bound 
states should be removed before the analysis is performed. 
Furthermore, the approximated double continuum states used 
here are not eigenstates of helium as the electron-electron 
correlation is ignored. The wave packet must therefore be 
propagated after the interaction until it reaches the Coulomb 
zone |j29D, where the electron-electron interaction is negligi- 
ble. 

To calculate the double ionization probability, we project 
the final wave function onto the non-correlated double contin- 
uum functions b, and sum over all contributions. 



-^double — ^ ^ I b^i^ ^7^2 ' ^(-^) 



(17) 



determine the total cross section for the nonsequential two- 
photon double ionization process H8n43ll . 
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where Iq is the pulse intensity. The finite duration of the pulse 
is accounted for by Tg//, which for a two-photon process 
reads lillii 



eff 





\Ht)] 


f 

-oo 


[ \ 
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For a sine-squared pulse envelope Te// = 35T/128, while 
for a Gaussian envelope T^f f^g = tq /2-y/7r/21n2. 



Having found the double ionization probability, we may then 



E. Accuracy and numerical convergence 



The reliability of numerically calculated quantities must be 
checked carefully, and this is usually performed by varying 
the relevant numerical parameters and studying the result- 
ing changes. It is possible, from certain physical considera- 
tions, to obtain a reasonable a priori estimate for the values 
of some parameters. In other cases, simple numerical calcula- 
tions may be performed to get such estimates. For example, in 
the present case the number of photons absorbed in the system 
determines the maximum value of the total angular momen- 
tum quantum number L that must be included in the basis, 
and also the required radial box size, from an estimate of the 
energy of the ejected electrons. The remaining radial basis 
parameters mainly determine the quality of the ground state 
(see Table and the maximum photoelectron energy sup- 
ported. By solving the related one-dimensional radial eigen- 
value problem (Coulombic potential, Z — 2), estimates for 
these parameters can be obtained. Specifically, estimating the 
density of states from 1/Ai?„, AEn being the energy separa- 
tion between state number n and ?i + 1 in the discretized basis, 
and comparing with the known value, we may determine to 
what extent the density of states is correctly represented in the 
box, in the energy region of interest. Finally, by sampling the 
parameter space in the vicinity of the values thus estimated, 
a good indication as to whether the results are converged is 
obtained. 

Time step size must also be considered, since our implicit 
time integrator incurs a local error of order At^. Incidentally, 
this error is one order of magnitude smaller than that of the 
exponential propagator it approximates, but the constants are 
typically different, and may depend on all the other parame- 
ters except t. In any case, we have used a default time step 
At = 0.01 a.u. Halving the time step to 0.005 a.u. produced 
changes in our calculated quantities of less than 0.1%. 
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Quantity 


Calculated value 


Reference value 




fa u ■) 


(a u ) 


ne \i D ) 


-jL.yyjj uu / 




He {2^S) 


-2.145 971 


-2.145 974 


He (135) 


-2.175 229 


-2.175 229 


H- {l^S) 


-0.527 735 


0.527 751 


Li+ {l^S) 


-7.279 827 


-7.279 913 



TABLE I: Calculated energy levels in the helium atom and the 
helium-like ions H~ and Li"'', compared with reference values from 
Drake [2]. The calculations used Imax = 7 and 50 exponentially 
distributed B-splines for each radial direction, in a box extending to 
50 a.u. 



III. RESULTS 
A. One photon double ionization 

Because of the electron-electron interaction, double ioniza- 
tion of helium may proceed through the absorption of a sin- 
gle photon. This one-photon double ionization process has 
been investigated at length in both theoretical 12. [H Htl and 
experimental [3, 113. l45ll studies, resulting in close quantita- 
tive agreement in both total and differential ionization cross 
sections, and a thorough understanding of the physical under- 
pinnings. This makes it an ideal benchmark against which 
new numerical schemes may be gauged. Accordingly, we 
have calculated cross sections for selected photon energies 
in the interval 80 — 200 eV. A box with Vmax — 160 a.u., 
311 B-splines, Imax = 5, and values of Lmax up to five was 
used. The pulse duration was set to 20 optical cycles, and the 
wave packet was propagated four additional cycles after the 
pulse was over, allowing the ionized component to reach the 
Coulomb zone. The results, shown in Fig. [T] include double 
ionization cross sections (red squares) and the ratio of double 
to single ionization (red circles), both within 5% of the exper- 
imental results of Samson et al. who state the accuracy of 
their results to be ±2%. 



B. Two photon double ionization 

We now consider the problem of two-photon direct double 
ionization. Our results, together with a small subset of re- 
sults from the numerous studies available in the literature, are 
shown in Fig. |2] The calculations have been made using a 
box extending to r^ax = 250 a.u., and 246 B-splines, while 
the angular basis was truncated at lmax — 5 and Lmax = 3, 
respectively. Additionally, the intensity of the laser field was 
fixed at 10^3 W/cm^, which is well within the perturbative 
regime. We also checked that decreasing the intensity by 
a factor of ten did not produce any significant changes in 
the results; at u = 51.7 eV the change in total cross sec- 
tion was less than 0.1%. Improving our basis by increasing 
{k,maxj2,max, Lmax) to (7,7,5) and the number of B-splines 
to 270 resulted in only minor changes in the cross sections 
(< 0.03% at 42 eV). To facilitate the separation of single- 
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FIG. 1 : (Color online) A comparison of calculated one-photon dou- 
ble ionization cross sections (red squares) and experimental values 
from Samson et al. if?] (blue full line). The ratio of double to sin- 
gle ionization is also shown. Red dots; calculated result, and green 
dashed line; experimental result. 

and double continuum components of the wave function, we 
ran the propagation algorithm an additional femtosecond af- 
ter the end of the pulse before performing the projections on 
the Coulomb waves, to ensure that the major part of the wave 
packet had entered the asymptotic Coulomb zone 1 18, 29]. 

The agreement between our results and those of Feist et 
al. [18] is particularly close, but this is not surprising due 
to the similarity of the numerical methods and the projection 
method used to extract the double ionization probability. In 
contrast, the J-matrix result of Foumouo et al. [5] (green line 
with diamonds in Fig.[2]i, and the perturbation theory result of 
Nikolopoulos et al. [17] (black line with triangles in Fig. |2j, 
deviate significantly from ours, i.e., the calculated total cross 
sections for the reaction differ by as much as an order of mag- 
nitude. In both of these approaches, correlation effects were 
included in the final state, to some extent, while in our calcula- 
tions no such effects were included. It should be noted, how- 
ever, that Foumouo et al. [5] obtained similar results to ours 
when they neglected completely the role of electron-electron 
interactions in the final wave function (black line with crosses 
in Fig.|2]i. This seems to stress the importance of electron cor- 
relations in the final states, however, by propagating the wave 
packet for a long time after the pulse, the correlation effect 
should, in prin ciple, be minimized, as argued and tested by 
Feist et al. 01 811 . In that particular study, they employed a very 
large grid and propagated the wave packet some 20 fs after 
the pulse, to explicitly check for convergence of the cross sec- 
tion. They also extracted the double ionization directly from 
the grid representation, by partitioning the radial grid, and by 
varying the partition limits found an upper bound for the pos- 
sible value of the cross section '--^25% higher than their quoted 
results (Fig.|2]i. 
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FIG. 2: (Color online) Two-photon double ionization cross sections. 
Blue line with circles: the present results obtained with Pyprop, black 
circle: experimental result of Hasegawa et al. 1 23], red cross: exper- 
imental result of Sorokin et al. |25], red line with squares: Feist et 
al. 1 18], green line with diamonds: Foumouo et al. |5] (with corre- 
lation, WC), black line with crosses: Foumouo et al. [5J (no corre- 
lation, NC), and black line with triangles: Nikolopoulos et al. [vh . 
The vertical lines define the two-photon direct double ionization re- 
gion. 



Regarding the apparent rapid rise in the value of the cal- 
culated cross section near the sequential ionization threshold 
(54.5 eV), this is usually attributed to the bandwidth of the 
pulses used in time-dependent methods and an unwanted in- 
clusion of the sequential process EE Him a. Thus, 

due to the finite spectral width of the pulse, the sequential 
process cannot be completely separated from the nonsequen- 
tial one, even below threshold. Now, since the relative im- 
portance of the sequential process scales as T^, as opposed to 
the T-dependence of the nonsequential, attempting to extract 
a cross section when both processes are present would result 
in a divergent behavior. This problem can be circumvented by 
simply increasing the pulse duration in order to lower its band- 
width. Following this procedure, one avoids significant con- 
tribution from the sequential component up to some finite dis- 
tance from the upper threshold, but at a certain point the over- 
lap with the sequential region again becomes non-negligible, 
and the scaling law causes an even sharper rise, due to the 
now longer pulse duration. Examining the relative importance 
of the different spectral components in the laser pulse as the 
pulse duration is increased, one can show that the relative con- 
tribution from the sequential process will ultimately become 
negligible, despite the T^-dependence of the ionization yield. 
Thus, using successively longer pulses, one may, at least in 
principle, resolve the behavior of the direct two-photon dou- 
ble ionization cross section arbitrarily close to the threshold, 
without contamination from the sequential process. 

Pursuing this line of thought, we have performed some ad- 



FIG. 3: (Color online) Fourier spectra of pulses with different tem- 
poral shapes and durations. Full lines represent sine-squared pulses 
with T = 2 fs (blue), 4 fs (black) and 6 fs (green), while the 
red dashed line represents a Gaussian pulse with to = 1.8 fs and 
T — 9.4 fs. The inset shows the 6 fs and Gaussian pulse spectra on 
a logarithmic scale. 



Pulse duration 
(fs) 


FWHM (fs) 


Cross section 


2 


0.7 


2.6 


4 


1.5 


2.1 


5 


1.8 


2.0 


6 


2.2 


2.0 


9.4* 


1.8 


2.0 



TABLE II: Double ionization cross sections at Tilu — 51.7 eV. 
('Gaussian pulse). 



ditional calculations at ui — 1.9 a.u. (51.7 eV) with longer 
pulse durations and different temporal shapes, whose Fourier 
spectra are shown in Fig. |3] The 2 fs pulse has a clear 
extension into the sequential region, indicated by the black 
vertical line, while the longer pulses have successively less 
overlap. The resulting cross section values are shown in Ta- 
ble |II] For pulse durations in the region 2 — 6 fs, and a sine- 
squared envelope, we found that the cross section leveled out 
at 2.0 X 10^^^ cm^s as T increased. Changing the temporal 
profile of the pulse to a Gaussian one, cf. Eq. |5] with a (to- 
tal) duration of ^ 9.4 fs, the same value for the cross section 
was obtained. In order to extract the cross section at photon 
energies exceeding 52 eV, significantly longer pulse durations 
than 5 fs would be required. Note that it was necessary to in- 
crease the radial box size to 350 a.u. (339 B-splines) in order 
to contain the double continuum wave packet when the pulse 
duration exceeded 4 fs. 
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IV. CONCLUSIONS 

In this paper we have presented a numerical method for 
solving the two-electron time-dependent Schrodinger equa- 
tion. After establishing the capability of the method through 
convergence tests and accurate reproduction of known phys- 
ical quantities, we applied the method to the study of two- 
photon direct double ionization of helium. Good agreement 
with several recently published results was found for the to- 
tal cross section. Finally, we investigated the behavior of the 
cross section near the sequential threshold, where a steep in- 
crease was observed. Calculating the value of the cross sec- 
tion at fixed frequency (5 1 .7 e V) for varying pulse duration, 
we found that the value converged for longer pulses, and it ap- 



pears that the cross section indeed exhibits a moderate grow- 
ing trend towards the threshold. 
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